Characterizing dynamics with covariant Lyapunov vectors 
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A general method to determine covariant Lyapunov vectors in both discrete- and continuous-time 
dynamical systems is introduced. This allows to address fundamental questions such as the degree 
of hyperbolicity, which can be quantified in terms of the transversality of these intrinsic vectors. For 
spatially extended systems, the covariant Lyapunov vectors have localization properties and spatial 
Fourier spectra qualitatively different from those composing the orthonormalized basis obtained in 
the standard procedure used to calculate the Lyapunov exponents. 
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Measuring Lyapunov exponents (LEs) is a central is- 
sue in the investigation of chaotic dynamical systems be- 
cause they are intrinsic observables that allow to quan- 
tify a number of different physical properties such as 
sensitivity to initial conditions, local entropy production 
and attractor dimension []]]. Moreover, in the context 
of spatiotemporal chaos, the very existence of a well- 
defined Lyapunov spectrum in the thermodynamic limit 
is a proof of the extensivity of chaos [2| , and it has been 
speculated that the small exponents contain information 
on the "hydrodynamic" modes of the dynamics (e.g., see 
[H and references therein). 

In this latter perspective, a growing interest has been 
devoted not only to the LEs but also to some corre- 
sponding vectors, with the motivation that they could 
contribute to identifying both the real-space structure 
of collective modes [4| and the regions characterized by 
stronger/ weaker instabilities [f|. However, the only avail- 
able approach so far is based on the vectors yielded by 
the standard procedure used to calculate the LEs p. 
This allows to identify the most expanding subspaces, but 
has the drawback that these vectors — that we shall call 
Gram-Schmidt vectors (GSV) after the procedure used — 
are, by construction, orthogonal, even where stable and 
unstable manifolds are nearly tangent. Moreover, GSV 
are not invariant under time reversal, and they are not 
covariant, i.e. the GSV at a given phase-space point are 
not mapped by the linearized dynamics into the GSV of 
the forward images of this point. 

While the existence, for invertible dynamics, of a 
coordinate-independent, local decomposition of phase 
space into covariant Lyapunov directions — the so-called 
Oseledec splitting [l[ — has been discussed by Ruelle long 
ago 0], it received almost no attention in the literature, 
because of the absence of algorithms to practically de- 
termine it. In this Letter, we propose an innovative ap- 
proach based on both forward and backward iterations 
of the tangent dynamics, which allows determining a set 
of directions at each point of phase space that are invari- 
ant under time reversal and covariant with the dynam- 



ics. We argue that, for any invertible dynamical system, 
the intrinsic tangent space decomposition introduced by 
these covariant Lyapunov vectors (CLV) coincides with 
the Oseledec splitting. 

As a first important and general application of the 
CLV, we show that they allow to quantify the degree 
of hyperbolicity of the dynamics. Considering that all 
physically relevant dynamical systems are not hyperbolic 
(i.e. stable and unstable manifolds are not everywhere 
transversal), and that many of the available theoretical 
results have been derived under the assumption of strict 
hyperbolicity (a prominent example being the Gallavotti- 
Cohen fluctuation theorem it is indeed highly desir- 
able to develop a tool to quantify deviations from hy- 
perbolicity. At the moment, this is doable only in very 
simple systems such as the Henon map or the Duffing os- 
cillator, where homoclinic tangencies can be detected by 
iterating separately the tangent dynamics forward and 
backward in time. Since CLV correspond to the local ex- 
panding/contracting directions, we can straightforwardly 
evaluate their relative transversality and, accordingly, 
quantify the degree of hyperbolicity. Note that GSV, 
being mutually orthogonal, are useless in this context. 
In a second important application of CLV we show that, 
contrary to the weak localization of GSV, they are gener- 
ically localized in physical space, providing an intrinsic, 
hierarchical decomposition of spatiotemporal chaos. Fur- 
thermore, the knowledge of CLV paves the way to analyt- 
ical methods for determining the LEs as ensemble- rather 
than time-averages. 

Description of the algorithm. We first summarize the 
standard method for computing the LEs (we consider, 
for simplicity, a V-dimensional discrete-time dynamical 
system). Let x n _i G 1Z N denote the phase-space point 
at time t n -\ and let {g^_x}, j = 1, • • • N, be the N or- 
thogonal vectors obtained by applying the Gram-Schmidt 
orthogonalization procedure to N tangent-space vectors 
(we shall call this the (n — l)th GS basis). Iterating 
the evolution equations once, g J n _ 1 is transformed into 
g J n = Jn-ign-i, where J„ is the Jacobian of the trans- 
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formation evaluated at time t n . The nth GS basis is 
thereby obtained by applying the Gram-Schmidt trans- 
formation to the vectors g^. This amounts to com- 
puting the so-called QR decomposition of the matrix 
G n = (g,\ | • ■ ■ |g^ ) whose columns are the Jacobian- 
iterated vectors of the (n — l)th GS basis: G„ = Q„R„. 
The nth GS basis is given by the columns of the orthog- 
onal matrix Q n , while R„ is an upper-triangular matrix 
whose off-diagonal nonzero elements are obtained by pro- 
jecting each vector g 3 n onto the subspace spanned by {g^j} 
with k < j. It has been shown 9] that, by repeating the 
above procedure up to a time t m for m much larger than 
n, the GS basis converges to an orthogonal set of vectors 
{e^}, k — 1, . . . , N - the mth Gram-Schmidt vectors - 
which solely depend on the phase space point x m . 

The LEs Ai > A2 > . . . > \n are then nothing but 
the time-averaged values of the logarithms of the diag- 
onal elements of R„. The method we propose also ex- 
ploits the usually disregarded information contained in 
the off-diagonal elements. Let us now assume that a set 
of GSV has been generated by iterating the generic ini- 
tial condition Xo. Let be a generic vector inside the 
subspace spanned by {e^}, k = l,...,j, i.e. the first 
j GSV at time t m . We now iterate this vector backward 
in time by inverting the upper-triangular matrix R m : if 
the = (e l m ■ u J m ) are the coefficients expressing it in 
terms of the GSV in x m , one has c^_ x = X) JRm]i/c lc m ; 
where [R]y is a matrix element of R. Since R m is upper- 
triangular, it is easy to verify that 6 S} t at all times t n . 
This is due to the fact that S 3 n is a covariant subspace. 
Iterating u J m backward for a sufficiently large number 
(rn — n) of times, it eventually aligns with the (back- 
ward) most expanding direction within S^. This defines 
v^, our intrinsic j-th (forward) expanding direction at 
the phase-space point x„. It is straightforward to verify 
that is covariant. Define the matrix [C m ]^ = c)° n \ 
then one has C m = R I „C m _ 1 . By multiplying both 
sides by Q m and substituting G m for its QR decompo- 
sition on the resulting right hand side, one is simply left 
with v?„ = J m _ 1 v^_ 1 for j = 1, . . . , N. The CLV are 
independent of where the backward evolution is started 
along a given trajectory, provided that it is sufficiently 
far in the future. Moreover, we have verified that they 
are invariant under time reversal, i.e. that the direction 
of v 3 n is the same whether we first move backward along 
a given trajectory (applying the standard orthonormal- 
ization procedure) and then forward (according to the 
above outlined methodology). 

Our CLV {v^j} thus constitute an intrinsic, covariant 
basis defining expanding/contracting directions in phase 
space [loj . The LEs are simply obtained from the CLV: 
the ith exponent is the average of the growth rate of the 
iih vector [11]. We have checked on simple invertible 
maps that they coincide with the Oseledec splitting in 
x m . We conjecture that this is the case for any invertible 
system. Note that our CLV are also well defined for non- 



invertible dynamics, since it is necessary and sufficient 
to follow backward a trajectory previously generated for- 
ward in time. In this respect they provide an extension 
of the Oseledec splitting. Finally, and retrospectively, a 
preliminary evidence of the validity of our approach was 
given in [12J, where CLV were introduced to character- 
ize time periodic orbits in a ID lattice of coupled maps. 
There, it was found that the number of nodes (changes 
of sign) in a CLV is directly connected to the position of 
the corresponding LE within the Lyapunov spectrum. 

We stress that the determination of the CLV can be 
very efficient, making them a truly practical tool (as op- 
posed, say, to calculating directly the Oseledec splitting 
in the case of invertible dynamics). Indeed, the major 
computational bottleneck is the memory required to store 
the matrices R n and the n-time GSV during the forward 
integration. This difficulty can be substantially reduced 
by occasionally storing the instantaneous configuration in 
real and tangent space and re-generating the rest when 
needed. 

Numerical analysis. We measured the CLV in four one- 
dimensional systems made of L nonlinear units coupled 
to their nearest neighbors. Periodic boundary conditions 
are used. The first is a chain of chaotic tent maps (TM) 
on the unit interval, 

4 +1 = (1 - 2e)/(<) + s [/(< +1 ) + /(xjr 1 )] 

with f(x) = ax if x < 1/a ^ 
and f(x) — a< "iZa^ otherwise. 

In the following we fix e = 0.2 and a = 2.3. 

The second system is a chain of symplectic maps (SM), 

Pn+1 = Pn + M ~ 4d ~ 9{<£n ~ ^rV 1 )] ft) 

Qn+1 = In + Pn+l 

where g(z) = sin(27rz)/(27r). This model was studied in 
[l3] ] to analyse the so-called "hydrodynamic Lyapunov 
modes". Eq. © conserves total momentum P = J2iP* > 
and is invariant under a translation of the q coordinates. 
Therefore, the Lyapunov spectrum possesses two null ex- 
ponents. In the following we fix fx = 0.6. 

The last two models are second-order continuous-time 
systems governed by 

qi = F(q i+ i-qi)-F(qi-qi-i). (3) 

For F(x) = sin(x), we have the rotator model (RM), 
while for F(x) — x + x 3 , the system reduces to a 
Fermi Pasta Ulam chain (FPU). These two widely stud- 
ied Hamiltonian systems provide a good testing ground to 
investigate the connection between microscopic dynamics 
and statistical mechanics. Besides the zero LE associated 
with a shift along the trajectory, both models have three 
other null LEs arising from energy and momentum con- 
servation plus translational invariance. Numerical simu- 
lations have been performed at energy density E/L = 1 
(for the RM) and E/L = 10 (for FPU). 



3 




<3> 7i/2 



7T/4 4> 



FIG. 1: (Color online). Probability distribution of the an- 
gle between stable and unstable manifold, (a) Henon map 
Xn+l = 1 — 1.4 x\ + 0. 3:Tn_i (green light line), and Lozi map 
x n +i = 1 — 1.4 \x n \ + 0.3x„_i (black line, rescaled by a factor 
10). (b) TM (L = 12, black dotted line), SM (L = 10, green 
dashed line), RM (L = 32, red dot-dashed line), and FPU 
(L = 32, blue full line). 



Hyperbolicity. A dynamical system is said to be hy- 
perbolic if its phase space has no homoclinic tangencies, 
i.e. the stable and unstable manifolds are everywhere 
transversal to each other. In the mathematical litera- 
ture, it is known that the Oseledec splitting is connected 
to hyperbolicity [3| , but the lack of practical algorithms 
to determine the splitting makes such results of little use 
in physically relevant contexts. Here, the knowledge of 
the CLV allows testing hyperbolicity by determining the 
angle between each pair (j, k) of expanding (J) and con- 
tracting (fc) directions 



■ 1 (|v^-vJt|)G[0,7r/2] 



(4) 



where the absolute value is taken because signs are irrel- 
evant. As a first test, we have computed the probability 
distribution P(</>) of 0*' for two classic two-dimensional 
maps. Arbitrarily small angles are found for the Henon 
map, while the distribution is bounded away from zero 
in the Lozi map (Fig. [TJi). This is perfectly consistent 
with the well-known fact that only the latter model is 
hyperbolic [lij ]. 

In spatially extended systems, given the multi- 
dimensional character of the invariant manifolds, it is 
appropriate to determine the minimum angle, <f>„ = 
minimi ( v « G E t , v « G E n)} wh ere £± are the ex- 
panding and contracting invariant subbundles at time t n 
along the trajectory. The histograms in Fig. [Tp show 
that models ([1]) and (|2|) are characterized by stronger 
hyperbolicity violations than the Hamiltonian systems. 
Altogether, recalling that $ refers to the least transver- 
sal pair of directions, we are led to conclude that the dy- 
namics of high-dimensional systems should be closer to 
hyperbolic than that of low-dimensional ones. This jus- 
tifies the often-made assumption that spatially-extended 
systems are practically hyperbolic. 

Localization properties in extended systems. The spa- 
tial structure of the vectors associated to the LEs is of 
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FIG. 2: (Color online). Inverse participation ratio Yz (see 
text) of both CLV and GSV for different dynamics. Time 
averages were performed over typically 10 5 ~ 10 6 timesteps 
and cubic splines have been employed to interpolate Y2 (h, N) 
between the discrete set of values h, j = 1, . . . , N. (a — c): 
Log-log plot of Y2 as a function of chain length L at fixed 
spectrum position h. CLV results are shown in full symbols, 
while GSV by empty symbols. In the log-log scale insets: 
inverse of the localization length £ has been subtracted from 
Y 2 to better show the CVL behavior Y 2 (L) ~l/£+ L~ 7 (see 
text). The dashed black lines mark a decay as 7 = i. (a): TM 
for ft = 0.1 (black circles) and h = 0.4 (red squares), (b): SM 
for h — 0.2 (black circles) and h = 0.4 (red squares), (c): FPU 
(h = 0.2, black circles) and RM (h = 0.2, red squares), (d): 
Lin-log plot of the asymptotic localization length £ of CLV as 
a function of h for TM (black circles) SM (red squares) and 
RM (blue triangles). 



interest in many contexts. We now show that the GSV 
— which have been used so far — and the CLV have quali- 
tatively different localization properties. One usually con- 
siders the inverse participation ratio [l6| ^2 = (X)i( a i) 4 ) 
where (•) indicates an average over the trajectory and a\ 
is a measure of the component of the jth vector at site i 
(with the normalization |a^| 2 = 1). In systems char- 
acterized by a single local real variable (such as our TM), 
a\ is taken to be the z-th component of the j-th CLV or 
GSV, while in the case of symplectic systems, where two 
components are present (v J = (Sq 3 , Sp 1 )), it is natural 
to choose («i) 2 = (5qf ) 2 + (Spj) 2 - In order to investigate 
the thermodynamic limit, it is necessary to determine 
Y2(h,L) for fixed h = (j — i)/L and increasing L. On the 
one hand, localized vectors are characterized by a finite 
inverse participation ratio, Yi{h^V) — > l/£, for L — > 00, 
where I is a localization "length" . On the other hand, in 
completely delocalized structures, Yi{h, L) ~ 1/L. 

In Fig. [5]we show how Y% typically scales with the chain 
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FIG. 3: (Color online). Trajectory averaged power spectrum 
(as a function of the wavenumber k = j 2tt/L, j = 1, . . . , L/2) 
of the space components of CLV (a) and GSV (b) correspond- 
ing to the smallest positive LE. Solid (black), dashed (red) 
and dot-dashed lines (blue) refer to FPU, RM and SM re- 
spectively (L — 512). The dotted green line, corresponding 
to a 1/k behavior is plotted for comparison in panel (a). 



and unstable manifolds are made available in generic 
models, many questions can be addressed in a more ac- 
curate way: Quantifying (non-)hyperbolicity in the con- 
text of the (numerical) attempts to "verify" the fluctu- 
ation theorem is one. Another set of questions relates 
to the spatial structure of the dynamics in extended sys- 
tems, such as the quantification of local degree of chaos 
(amount of instability), a hierarchical decomposition of 
spatiotemporal chaos, the search for true, intrinsic, col- 
lective ( "hydrodynamic" ) modes, etc. A further field 
where the knowledge of CLV can help to make progress 
is optimal forecast in nonlinear models. Here the knowl- 
edge of the local transversality of the invariant manifolds 
can indeed be combined with the so-called bred vectors 
to use the information on the past evolution to decrease 
the uncertainty along unstable directions [l8| . 



length L. The GSV show weak (de)localization: their 
participation ratio exibits an /i-dependent "dimension" 
"q(h): Yi ~ L~ n ( h \ One can show that this anomalous 
behavior is entirely due to the Gram-Schmidt procedure, 



and has nothing to do with the dynamics [17J. On the 
other hand, CLV are localized objects. For TM, SM and 
RM dynamics we find good evidence of the scaling law 
Y 2 (h, L) ~ l/l{h)+L-"i with 7«i. This allows for a re- 
liable determination of I. For the FPU dynamics, we find 
only slight curvature in the log-log plot of Fig. sig- 
nalling that larger system sizes are probably needed to 
definitely enter the scaling regime. Moreover, for sym- 
plectic dynamics the localization length 1(h) diverges as 
h — ► 1 (Fig. Assuming the continuity of the LE 

spectrum, the divergence of I is not surprising, since the 
conservation laws imply that the Lyapunov vectors (both 
GSV and CLV) corresponding to h — 1 (i.e. to null LEs) 
are completely delocalized. 

Fourier analysis. Another way proposed to character- 
ize the spatial structure of a Lyapunov vector is to look 
at its power spectrum S(k) = | J2 m An eimfe | > where /3 m 
denotes the vector component associated with the space 
coordinate q m at site m. For instance, this was used in 
[l3| in the context of the investigation of so-called "hy- 
drodynamic" modes (only GSV were considered there). 
Here, we have focused on the vector corresponding to the 
smallest positive LE in our symplectic models, for which 
this LE goes continuously to zero as the system size in- 
creases (note that GSV and CLV coincide for the null 
exponents linked to symmetries and conservation laws). 
We observe again a clear qualitative difference between 
the spectra of GSV and CLV (Fig. [3]). In particular, the 
near-zero CLV exhibit an intriguing low-frequency diver- 
gence of the 1/k type in all three symplectic models we 
have analysed. Thus, the qualitative difference between 
GSV and CLV extends to the h — > 1 case. 

Perspectives. Now that the local directions of stable 
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